function [rho,rhou,E] = MacCormack3(n,v)
%用LaxFridrichs格式计算Euler方程的初边值问题3
%参数说明：
%   输入:
%   n：单位长度上的网格数量。
%   v:\tau / h
%   输出：
%   rho,rhou,E:相应变量在时刻T的值。
    T = 0.381;
    h = 1/n;
    tau = h*v;
    step = floor(T/tau);
    %初始化
    rho = zeros(n+1,2);
    rhou = zeros(n+1,2);
    E = zeros(n+1,2);
    
    %初始条件
    rho(:,1)=1;
    rhou(:,1)=0;
    a = floor(0.1*n);
    b = floor(0.9*n);
    %    E(:,1)=sin(4*pi/(n+1)*[1:n+1]);
    E(1:a,1)=2.5*10^3;
    E(a+1:b,1)=2.5*10^-2;
    E(b+1:n+1,1)=2.5*10^2;

    %边界条件

    rho(1,2)=rho(1,1);
    rho(n+1,2)=rho(n+1,1);
    E(1,2)=E(1,1);
    E(n+1,2)=E(n+1,1);

    
    frho=zeros(n+1,1);
    frhou=zeros(n+1,1);
    fE=zeros(n+1,1);
    
    %定义存储中间步的变量。
    irho=zeros(n+1,1);
    irhou=zeros(n+1,1);
    iE=zeros(n+1,1);
    
    irho(1)=rho(1,1);
    irho(n+1)=rho(n+1,1);
    irhou(1)=0;
    irhou(n+1)=0;
    iE(1)=E(1,1);
    iE(n+1)=E(n+1,1);

%    plot(E(:,1));
%    
%    axis([0,size(E,1),min(E(:,1)),max(E(:,1))])
%    pause()
%

    %迭代
    ax = 0:1/n:1;
    ti = 2;
    t = 1;
    for s = 2:step+1
        x = ti;
        ti = t;
        t = x;
        %计算f_j

        rhou2 =rhou(1:n+1,ti).^2./rho(1:n+1,ti);%rho*u^2
        p = 0.4*(E(1:n+1,ti)-rhou2/2);
        
        frho(1:n+1)=rhou(1:n+1,ti);
        frhou(1:n+1)=rhou2+p;
        fE(1:n+1)=rhou(1:n+1,ti)./rho(1:n+1,ti).*(E(1:n+1,ti)+p);

        %计算中间步

        irho(2:n)=rho(2:n,ti) - v*(frho(3:n+1)-frho(2:n));
        irhou(2:n)=rhou(2:n,ti) - v*(frhou(3:n+1)-frhou(2:n));
        iE(2:n)=E(2:n,ti)  - v*(fE(3:n+1)-fE(2:n));

        % 计算中间步的f_j
        for i = 1:n+1
            rhou2 =irhou(i)^2/irho(i);
            p = 0.4*(iE(i)-rhou2/2);
            frho(i)=irhou(i);
            frhou(i)=rhou2+p;
            fE(i)=irhou(i)/irho(i)*(iE(i)+p);
        end        
        %下一步
        for i=2:n
            rho(i,t) = (( rho(i,ti)+irho(i) ) - v*(frho(i)-frho(i-1)))/2;
            rhou(i,t) = (( rhou(i,ti)+irhou(i) ) - v*(frhou(i)-frhou(i-1)))/2;
            E(i,t)= (( E(i,ti)+iE(i) ) - v*(fE(i)-fE(i-1)))/2;
        end   
        
        plot(ax,rho(:,t));
        pause()
    end        
end